Thermodynamic Entropy and Chaos in a Discrete Hydrodynamical System 



Franco Bagnoli* 
Dipartimento di Energetica, Universita di Firenze, 
Via S. Marta 3, 1-50139 Firenze, Italy 
also CSDC and INFN, sez. Firenze 

Raul Rechtmant 
Centro de Investigation en Energia, Universidad National Autonoma 
de Mexico, Apdo. Postal 34, 62580 Temixco, Mor., Mexico 
(Dated: January 10, 2009) 

We show that the thermodynamic entropy density is proportional to the largest Lyapunov ex- 
ponent (LLE) of a discrete hydrodynamical system, a deterministic two-dimensional lattice gas 
automaton. The definition of the LLE for cellular automata is based on the concept of Boolean 
derivatives and is formally equivalent to that of continuous dynamical systems. This relation is jus- 
tified using a Markovian model. In an irreversible process with an initial density difference between 
both halves of the system, we find that Boltzmann's H function is linearly related to the expansion 
factor of the LLE, although the latter is more sensitive to the presence of traveling waves. 
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I. INTRODUCTION 

The relation between thermodynamics and the under- 
lying chaotic properties of a system is of great relevance 
in the foundations of statistical mechanics [1, 2] and has 
attracted much interest. What is the relation between 
chaotic dynamics and thermodynamics? Is chaos rele- 
vant for thermodynamics? Some interesting results have 
been found for some models. In the case of a family of 
simple liquids, including Lennard- Jones, a relation has 
been found between the Kolmogorov-Sinai entropy and 
the thermodynamic entropy [3] . Lorentz gases have been 
extensively studied and relationships between chaotic dy- 
namical properties and transport coefficients have been 
established [4-8]. The discussion of this problem for 
other simple models and maps is extensive [9-13]. 

Lorentz gases are essentially one-particle systems. In 
this paper we present a simple model of a gas with a large 
number of interacting particles and find interesting re- 
lations between dynamic and thermodynamic quantities, 
somewhat motivated by Bernoulli systems. For these, the 
Kolmogorov-Sinai entropy is on the one hand the sum of 
the positive Lyapunov exponents, and on the other, the 
thermodynamic entropy scaled by a time constant corre- 
sponding to the correlation time [14, 15]. 

The model we study is a deterministic lattice gas 
cellular automaton (LGCA). LGCA are simple models 
with hydrodynamical behavior [16-18]. In particular, the 
D2Q9 LGCA is a two dimensional model with nine veloc- 
ities and is one of the simplest models where equilibrium 
thermodynamics can be found with a non trivial temper- 
ature [19]. For cellular automata, Lyapunov exponents 
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can be defined using Boolean derivatives [20] . This defi- 
nition is formally similar to that of continuous dynamical 
systems, a long time average of the logarithm of the lin- 
earized expansion factor [21]. We also look for a relation 
between Boltzmann's H function in a simple irreversible 
process, and the logarithm of the expansion factor. 

The paper is organized as follows. In Sec. II we present 
a deterministic version of the D2Q9 model together with 
an introduction to Lyapunov exponents of cellular au- 
tomata. In Sec. Ill we discuss the equilibrium thermody- 
namics of the D2Q9 models and in Sec. IV we show there 
is a close relationship between the equilibrium entropy of 
the model and its largest Lyapunov exponent. This re- 
lationship is explained by finding the Kolmogorov-Sinai 
entropy of a Markov chain which relates the thermody- 
namic entropy to the largest Lyapunov exponents. We 
also show that Boltzmann's H function goes to its equi- 
librium value in an irreversible process in the same way 
the Lyapunov exponent does. We end with a discussion 
on why these quantities are related. 



II. THE D2Q9 MODEL 

The D2Q9 model is defined on a two dimensional 
square lattice [22]. The evolution is in discrete time 
steps where unit mass particles at every site r can move 
with one of nine velocities Cq = (0,0), C\ ~ (1,0), 
c 2 = (0,1), c 3 = (-1,0), c 4 = (0,-1), c 5 = (1,1), 
c 6 = (-1,1), c 7 = (-1,-1), c 8 = (1,-1). The state 
of the automaton is given by the set of occupation num- 
bers s(t) = {sk{r,t)} : where Sk(r,t) — 1 (0) indicates 
the presence (absence) of a particle with velocity at 
site r and time t. An exclusion principle forbids the pres- 
ence of more than one particle in a given site and a given 
time with a given velocity. 

The time evolution of the system is given by collision 



FIG. 1: All these states have the same number of particles, 
momentum and energy. An arrow represents the presence of 
a particle with the velocity of the arrow, the open circle a 
particle at rest. 



and streaming operations. In the collision step, the parti- 
cles at every site change their velocities conserving mass, 
momentum and energy. In the streaming operation par- 
ticles move to neighboring sites according to their veloc- 
ities. 

Since the number of states for a given site is finite (2 9 ), 
the local collision operator C is generally implemented 
as a look-up table. Given a local configuration, the con- 
servation constraints may not define the outgoing state 
completely as the example in Fig. 1 shows. If any one of 
the six states shown is chosen as the input state, any one 
of the other six states can be the output state. Therefore, 
the C look-up table has several columns for all the possi- 
ble output states. In order to make the automaton deter- 
ministic, we first find for all input states, the number of 
output states. We then find the least common multiplier 
of the number of output states and construct a new C 
look-up table with this number of columns. These are 
filled by repeating the output states the required num- 
ber of times. The least common multiplier is sixty, so 
the rows which have any one of the input states of Fig. 1 
are repeated ten times. At the beginning we assign an 
integer random number r](r) to each site, to be used in 
choosing the column from which the output state is cho- 
sen (quenched disorder). That is 

s k (r, t + 1) = C k {s {r, t),...,s 8 (r, t), V {r)). (1) 

The choice of the quenched disorder is analogous to the 
random disposition of scatterers in a wind-tree or other 
similar Lorentz gases [23, 24]. 

The evolution of the model can also be made reversible. 
In order to do so, the collision table must satisfy the 
condition that in every column ,s = ICIC (s) holds where 
I denotes the operator that inverts the velocities of any 
state s. This means that a collision C followed by an 
inversion I and another collision and inversion leaves the 
state unchanged when it is taken from the same column 
in the collision table as we show in Fig. 2. Once we have 
a look-up table that is deterministic as described above, 
the elements in each row are rearranged to satisfy the 
reversibility condition. 

In order to introduce the concept of Lyapunov expo- 
nents for such a system, it is convenient to simplify the 
notation used. We use the index n to indicate the posi- 
tion r and the velocity k, with n — 1, . . . , 9L, where L 
is the number of sites of the automaton. A configuration 
of the LGCA at a given time is given by 9L occupation 
numbers (bits), and its evolution can be seen as an ap- 
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FIG. 2: In this example, the input state at some fixed site is 
the one shown on the left of Fig. 1 and in the column number 
assigned to this site, the output state is the fourth one of that 
Fig. These are shown in the upper part of this Fig. Then, by 
inverting the velocities we get the state on the lower left and 
the collision table should contain in the same column for this 
input state the output state shown on the lower right that 
upon inversion of the velocities yields the original incoming 
state. 



plication of a set of Boolean functions 

s n (t + l) = F n (s(t)). (2) 

The functions F n represent different entries of the colli- 
sion table, and differ in the velocity index k and quenched 
disorder r}{r). However, since the distribution of the dis- 
order is uniform, and, as shown in the following, the cor- 
relations among variables decay very fast, the system is 
translationally invariant at a mesoscopic level. 

Let s(0) and x(0) be two initially close configurations, 
for example all the components of x(0) may be equal to 
those of s(0) except for one. We define the bitwise dif- 
ference between these two configurations with the term 
"damage" . The smallest possible damage is one and the 
damage vector v(0) has a one in the component where 
s(0) and x(0) differ and has zeroes in all the others. If 
this damage grows on average during time, the trajec- 
tory is unstable with respect to the smallest perturba- 
tion. However, due to the discrete nature of LGCA, de- 
fects may annihilate during time evolution, altering the 
measure of instability of trajectories. The correct way of 
testing for instability is that of considering all possible 
ways of inserting the smallest damage in a configuration, 
using as many replicas as the number of components of 
the configuration, and letting them evolve for one time 
step counting if the number of damages has grown or 
diminished. The ensemble of all possible replicas with 
one damage each is the equivalent of the tangent space 
of discrete systems. 

The task of computing the evolution in tangent space 
is clearly daunting, but by exploiting the concept of 
Boolean derivatives [20, 25] it is possible to develop a 
formula very similar to the one used in continuous sys- 
tems. The Boolean derivative is defined by 

= \F n (...,s p ,...)-F n (,..,l-s p ...)\ (3) 
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with n,p = 1, . . . , 9L. This quantity measures the sen- 
sitivity of the function F n with respect to a change in 
s p . The Jacobian matrix J(s) has components J np = 
dF n (s)/ds p . 

We now consider the map 



v(t + l) = J(s(t))v(t) 



(4) 



with v(0) as mentioned above. It is easy to check that 
\v(t)\ = J2i v i(t) 1S the number of different paths along 
which a damage may grow in tangent space during time 
evolution, i.e., with the prescription of just one defect per 
replica [20]. If there is sensitivity with respect to initial 
conditions, one expects that |v(T)|/|r>(0)| ~ cxp(A T T) 
for large T where A T is the largest finite time Lyapunov 
exponent (LLE). It then follows that 



T-l 



A t = 4 E 1o S u W = ( 1o S m )t 



(5) 



t=i 



with u(t) = \v(t)\/\v(t — 1)| the expansion factor of 
the LLE. The definition should include the limit when 
T — > oo but in practice we always evaluate the finite time 
LLE with T. The LLE depends in principle on the initial 
configuration s(0), initial damage v(0) and quenched dis- 
order 77, but in practice it assumes the same value for all 
trajectories corresponding to the same macroscopic ob- 
servables when T is sufficiently large. The LLE of CA as 
defined above has been used to classify elementary and 
totalistic Boolean cellular automata [20, 26]. 



III. THERMODYNAMICS OF THE D2Q9 
MODEL 



The single particle velocity distribution functions 
fk(r,t) are defined as the average number of particles 
at site r with velocity c k at time t over R samples that 
share the same quenched disorder rj and macroscopic con- 
straints with different microscopic initial configurations. 

We checked that the single site two-particle correlation 
function factorizes into the product of single particle dis- 
tributions before and after the collision, in equilibrium 
and out of equilibrium conditions. In equilibrium, the 
correlation function decays to zero for just one lattice 
spacing. Out of equilibrium, starting with very different 
configurations in the two halves of the system, the cor- 
relation function, still being quite small, exhibits a cor- 
relation length of some lattice spacings for a short time. 
This corresponds to the coherent motion of particles in 
a shock wave, where the local density is near to zero or 
nine. However, this correlation quickly disappears; al- 
though the motion is correlated at a macroscopic level, 
as soon as the local density of particles is different from 
the extremes (for which the collision table has few output 
configurations) the velocities quickly decorrelate. 

The thermodynamic entropy density can be found an- 



alytically in the thermodynamic limit as follows [27]. Let 
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with N and E the number of particles and total energy 
respectively, N k the number of particles in direction k, 
and eo = 0, £1,2,3,4 = 1/2, £5,6.7,8 = 1- The number 
density n, the equilibrium density functions f k and the 
energy density e are 



N 



n 



fk 



E 



(7) 
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with L the number of sites in the lattice. The micro 
canonical partition function Q is 



n = Jl[dN k S(N -Y,N k )S(E -J2^N k ): 



u(N ,...,N s ) 
with 5 the Dirac-i5 function and 



n 



N k \{L-N k )\ 



(8) 



(9) 



the number of microscopic states in Fermi-Dirac statis- 
tics. 

The thermodynamic entropy is S{E 1 N 1 L) = 
logQ(E,N,L). Using Stirling's approximation and the 
Fourier transform of the S function 



S = log L £ 



Y[df k J dp 1 dp 2 exp(La) 



(10) 



with 



= - E [A l0 § /* + (!" fk) lo g(! - fk)] 



+ ipi 



W2 



(11) 



The integrals of Eq, (10) can be evaluated using the sad- 
dle point method and the result is exact in the thermody- 
namic limit. The thermodynamic limit entropy density 
s is 

s(e,n) = lim —S(eL,nL,L) 

= ~ E lo § fk + (1 - h) l°g(l - h)] ■ (I 2 ) 

k 

In the last expression /o, . . . fg,pi, and pi maximize a 
given by Eq. (11). 

The quantities p\ and pi are related to temperature 
and chemical potential. Taking the partial derivative of 
a with respect to f k and equaling to zero we find that 



log 



fk 



1 - fk , 



= -ipi - ipid kl 



(13) 
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and that 

fk = [1 + cxp(zpi + ip 2 ek)r 1 ■ 
Using Eq. (13) the entropy density is 

* = -£(lGg(l-/fc) 



The Eulcr equation is 



_ e P iiN 

^ rp rp rp 



(14) 



(15) 



(16) 



where P is the pressure and \x is the chemical potential. 
Comparing the last two expressions 



A* 



= -5> g (i-/ fc ). 



(17) 
(18) 

(19) 



FIG. 3: Distribution functions in equilibrium as a function of 
e for n = 4. The continuous curves are the calculated values, 
the symbols are the results of numerical simulations. For 
e = ep there is a rest particle in all sites and no fast particles 
moving along the diagonals. The other three particles can 
occupy one of the four slow directions. For e = eu there are 
no rest or slow particles, the four particles are moving along 
the diagonals. 



Now, the equilibrium distribution functions are 
fk = [l + exp(e fe /T- At /T)]- 1 . 



(20) 



The distribution functions /o, /i, and f$ dehne the 
equilibrium state since fx = fa = = fi and fa — 
f% = fj = fs- These satisfy the conservation equations 
of mass end energy and maximize a. We also note that 
these distributions satisfy another condition equivalent 
to the latter. From Eq. (13) for k = 1,5 and for k = 
and Eq. (18) 



fo 



log 



fl 



1-/6 



Then 



/i 2 (l-/o-/ 5 )=/o/ 5 (l-2/ 1 ). 



(21) 



This last expression, together with mass and energy con- 
servation determine the equilibrium values of the distri- 
bution functions. In Fig. 3 we show the values of the 
equilibrium distribution functions as continuous curves, 
together with result of numerical simulations. The agree- 
ment between numerical simulations and computed val- 
ues of the probability distributions is perfect, and this 
confirms that the system is a perfect gas. Further inves- 
tigations on correlations in out of equilibrium conditions 
are reported in the following Section. 



IV. RESULTS AND DISCUSSION 

In Fig. 4 we show the entropy density s = S/L and the 
LLE A as a function of e for fixed n both calculated nu- 
merically. The entropy density is found by substituting 
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FIG. 4: Entropy density s (solid line) and LLE A (dashed 
line) for n = 4, simulations with R = 40 in a 40 x 40 lattice. 



the numerical values of fk, for example those of Fig. 3, in 
Eq. (12). The entropy density grows for ep < e < e* and 
then decreases for e* < e < cm with e* = 2n/3. The en- 
ergy density e* is the value for which / = /i = f$. The 
largest Lyapunov exponent A shows the same behavior 
having a maximum also at e* . The results shown suggest 
that s is proportional to A. This is emphasized in Fig. 5 
where the data lie near a straight line for different val- 
ues of n. This result shows that the LLE grows with the 
number of available states measured by the equilibrium 
entropy. 

The proportionality between the thermodynamic en- 
tropy density and the largest Lyapunov exponent can be 
understood in the framework of the stochastic approxi- 
mation of a chaotic dynamics [28] . In this approximation, 
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FIG. 5: We show s versus A for n — 2 (diamonds), n = 4 
(solid line), n = 4.5 (dotted line) and n — 6 (pluses). The 
simulations are performed for different values of e from eF (n) 
to eni(n). 

the dynamics is considered at discrete time intervals, and 
generates a discrete and finite partition of coarse-grained 
states labeled by an index i, with i = 1, . . . , M. The 
evolution is represented by a Markov chain where Wij 
is the transition probability from state j to state i and 
Si W»j = 1 . The asymptotic probability distribution 
is denoted by pj, with = 1 an( ^ 12j^ijPj = P»- 

The Kolmogorov-Sinai entropy K per time interval r 
is [21, 28] 

3 i 

and the entropy density s is 

s = -^^logpi. (23) 

i 

For Bernoulli systems the equilibrium is reached in just 
one time step, Wij = pj, and trivially K = t~ 1 s [29]. For 
a somewhat more general process, we assume that pi = 
l/M (microcanonical distribution) so that s = logM. 
We further assume that the transition matrix is irre- 
ducible (so that the system is ergodic) and that in each 
row and column there are ct^M non zero equal entries 
with values 1/ajM, < on < 1. Then 

K= (log a) + log M, (24) 

where (log a) is the average of logaj. 

If we further assume the validity of the Pesin rela- 
tion [30], K equals the sum of positive Lyapunov ex- 
ponents. For many Hamiltonian systems and symplcc- 
tic maps, the Lyapunov spectrum takes the form Aj ~ 
Ao(l —i/(2N)), where Aj is the i-th Lyapunov exponent, 
2 = 0,..., 2N — 1 and the number of particle iV is large. 
In these cases, the Kolmogorov-Sinai entropy per degree 
of freedom is proportional to the largest Lyapunov ex- 
ponent Aq. The proportionality constant, however, may 



FIG. 6: Boltzmann function H (thick line) and Lyapunov 
expansion factor (£} (oscillating thin line) versus time t, sim- 
ulations with til = 7.2, e_L = 4.8, tir — 1.8, e_R = 1.2, and 
R = 40 on a 40x40 lattice. The inset shows that, disregarding 
the oscillations of {£) , there is a linear relation between these 
quantities. The dashed line is the best fit H — 1.15{£) + 0.2. 

depend on the value of control parameters, namely the 
energy. 

The shape of the Lyapunov spectrum is roughly lin- 
ear for the product of random matrices with the struc- 
ture of (weakly) locally coupled Hamiltonian chaotic sys- 
tems [31] and for coupled nonlinear oscillators [32]. If we 
assume that also in other cases the shape of the spectrum 
(which in general is not linear, see for instance Ref. [33] 
for the hard spheres case) does not change with the en- 
ergy, it follows that s and A are linearly related. This 
last assumption is probably the weakest one, at least for 
LGCA for which the Lyapunov spectrum is ill-defined, 
and we think it is the main reason for the discrepancy 
from linearity in Figure 5. 

In our model, the value of the LLE is related to the 
number of ones in the Jacobian matrix J defined by 
Eq. (3). This Jacobian matrix contains the linearized 
effects of the streaming and collision operators. Stream- 
ing corresponds to a scrambling of the components of 
the tangent vector v(t) and therefore does not alter its 
norm. This is left to collisions, when more than one out- 
put configurations are possible. The number of equiva- 
lent output configurations in the collision table is small 
for the extreme values of number and energy densities, 
and larger for intermediate values. Similar considerations 
apply to the number of equivalent configurations for a 
given macroscopic distribution of density and velocities, 
and constitute the microscopic origin of the relation be- 
tween statistical and dynamical quantities. 

We now discuss an irreversible process where a square 
lattice is initially in an equilibrium state with the left and 
right sides having different number and energy densities 
til, ur, eL, and e^. The system evolves toward equilib- 
rium by means of damped traveling waves. The single 
site two-particle correlation function, although small, ex- 



6 



hibits a correlation length of some lattice spacings for a 
short time. 

Macroscopically, one observes a coherent motion of 
particles in a shock wave, where the local density is near 
zero or nine. Although the motion is correlated at a 
macroscopic level, as soon as the local density of parti- 
cles is different from the extremes (for which the collision 
table has few output configurations) the velocities quickly 
decorrclate. 

Boltzmann's H function is defined by 

H(t) = -J2fk(r,t)logf k (r,t). (25) 

r,k 

In the numerical simulations the distribution functions 
are averaged over R samples and the average Lyapunov 
expansion factor is (£) = (1/i?) X^iLi logu^L As we show 
in Fig. 6, the two quantities exhibit similar behavior. 
The Lyapunov expansion factor exhibits more marked 
oscillations, indicating that it is more sensible to the local 
variations in density. The inset of Fig. 6 shows that, 
disregarding oscillations, {£) is linearly related to H for 
all the relaxation phase. 

We can identify several time scales in this irreversible 
process. There is a fundamental time scale, which is fixed 
to unity. There is also a correlation time scale, which is 
of the order of the mean free time which depends on the 
occurrence of non-trivial collisions and is larger where 
the density is either small or large. During shock waves, 
locally one may have variations of the density and there- 
fore correlations, as already reported, of the order of some 
time steps. A third time scale is given by the oscillations 
induced by the traveling shock waves. This is a dynam- 
ical, mesoscopic scale, that depends on the size of the 
system, and is revealed by the oscillations of (£) shown 
in Fig. 6. The slowest time scale is given by the relaxation 
to equilibrium, shown both by H and by (£) . 



The computation of (£) is performed using a set of tan- 
gent vectors, Eq. (4), and these vectors constitute a sort 
of local memory of the past state. In systems with local 
variations of density, as in our system in the presence of 
traveling waves, statistical quantities like H depend on 
the instantaneous state of the system, while dynamical 
ones like {£) depend also on the variations of this state. 
This factor may be the origin of the different relation be- 
tween statistical and dynamical quantities in equilibrium 
and during the relaxation phase. 



V. FINAL REMARKS 

The D2Q9 reversible LGCA model we have discussed 
exhibits hydrodynamical and thermodynamical behavior. 
For this model, the thermodynamic entropy density is 
proportional to the largest Lyapunov exponent. Also, in 
a simple irreversible process, Boltzmann's H function is 
proportional to the expansion factor of the largest Lya- 
punov exponent. A simple stochastic coarse grained dy- 
namics can explain the proportionality between s and A. 
We note that this result does not depend on the D2Q9 
dynamics and should hold for other, possibly more realis- 
tic models. Finally, the fact that the thermodynamic en- 
tropy density and the LLE of the model are proportional 
confirms that the definition of the latter is appropriate. 
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